Finite size effects in the specific heat of glass-formers 
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Abstract. We report clear finite size effects in the specific heat and in the relaxation times of a model glass former at 
temperatures considerably smaller than the Mode Coupling transition. A crucial ingredient to reach this result is a new Monte 
Carlo algorithm which allows us to reduce the relaxation time by two order of magnitudes. These effects signal the existence 
of a large correlation length in static quantities. 



INTRODUCTION 

When compared to experiments, numeric simulations 
of glass-forming systems present severe limitations. In 
fact, while the actual physical origin of the phenomenon 
termed generically as glass transition is still a matter of 
debate, its investigation is limited by the practical impos- 
sibility to reach thermodynamic equilibrium when the 
relaxation time T of a sample glass in a laboratory be- 
comes of order of 10 2 seconds or when the T of a model 
glass in a modern computer becomes of ~ 10~ 8 seconds. 
Thus the glass temperature T g is conventionally defined 
in experiments and numeric simulations cannot get very 
close to it 1 . Furthermore, the statistical accuracy in ex- 
periments is fairly larger, because the time-length Jz? of 
typical experiments is many decades larger than in sim- 
ulations. Indeed, statistical errors decrease roughly as 
1 / V^ind where the number of independent configura- 
tions jVind is given by J/m ~ ££ ji. This problem is par- 
ticularly serious in the study of time correlators, whose 
peculiar behavior is one of the most clear signature of the 
glass transition. For every observable, O, the normalized 
time correlator is 
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The region C ~ is very interesting, since it allows to 
estimate the relaxation time (see below). On the other 
hand, in that region the relative error from a single mea- 
surement grows as C hence a very large JVi n d is needed 
to get some significant signal over the noise. As a conse- 
quence, most of previous numeric work was confined to 



e.g. in S1O2 Molecular Dynamics simulation have been performed 
only down to T • 



C > 0.1, while in experiments one is able to explore far 
lesser values of C. 

Despite those limitations, numeric simulations turn 
out to be quite useful in many respects, mostly in the 
key issue of identifying the quantity suffering the largest 
spatial fluctuations at the glass transition. Ijjyl |3l|4j|5|| 
In fact, unlike standard phase transitions the equilibra- 
tion time of glass-formers (supercooled liquids, poly- 
mers, proteins, superconductors, etc.) diverges without 
dramatic changes in their structural properties 1 6, 0, 0], 
challenging the interpretation of the slowing down as a 
critical phenomenon|9|. A possible way out has been 
found in the heterogeneous dynamic scenario, which 
postulate that local time correlators (i.e where O is a lo- 
cal quantity) are correlated over diverging length scales. 
The possible advantages of numeric simulations with re- 
spect of experiment are twofold. First, they are more suit- 
able to investigate local quantities. Second, they allow 
an easier implementation of Finite Size Scaling tech- 
niques, apowerful tool to detect diverging correlation 
length. I lOl fill Il2ll New, fast algorithms for numerical 
simulation are therefore needed in order to to get as close 
as possible to the experimental regime, while retaining 
all the advantages provided by numerics. 



TRANSLATIONAL INVARIANT 
QUANTITIES 

In this numeric work, we wish to ask the following 
question: Are there static observables whose correlation 
length 4 diverge at the glass transition? (Note that exper- 
imental clues are not available). In other words we won- 
der if a standard second order phase transition is behind 
the glass transition. We borrow from the dynamic het- 
erogeneity scenario the hypothesis that a diverging cor- 



relation length would show up in four-particle correlation 
function and focus on translational invariant quantities. 

Little attention has been paid so far to such quantities. 
In fact experiments (e.g. elastic scattering for the spatial 
correlations and dielectric response for the dynamic be- 
havior) and simulations typically investigate quantities 
which are not invariant under an uniform displacement 
of the coordinates, F, = + A, namely the density fluctu- 
ations of generic momentum Q: 

p(3) = ^£M ( 2 ) 

where p(0) is the particle density and it is constant 
because of mass conservation. Density fluctuations at 
small Q are surely slow modes of the system, their 
time-correlator decaying exponentially with characteris- 
tic time T ~ D l Q 2 (D is the diffusivity). We have no 
guarantee whatsoever though that they are the only rele- 
vant slow modes. 

How to expand the space of observables investigated 
taking into account even quantities which preserve the 
translational symmetry? For the sake of definiteness we 
limit to the NVT ensemble and introduce energy density 
fluctuations: 

where is the interaction energy V {f^ — ri) between 
particles k and j. Note that p e (0) is the non-conserved 
potential energy density, which we call e (the internal en- 
ergy is then \k-&T + (e)). A general way to obtain trans- 
lationally invariant observables consists in multiplying 
densities with wave-vectors q and — q: 

J?(q)=p(q)p(-q), S" e tf) = Pe(q)pe(-q) ■ (4) 

We see then that the correlation functions of these 
translational invariant quantities involve the positions 
of four particles. Moreover the mean value of S^{q) 
is the static structure factor S(q), while (^(0)) is re- 
lated to the constant-volume specific heat, Cy as T 2 Cy = 
N(o Q ) 6 [{y e (0))-(e) 2 } . 

THE LOCAL SWAP MONTE CARLO 
ALGORITHM 

We study a 50% mixture of particles interacting through 
the pair potential V a p{r) = e[(a a + Op)/r] n +C a p, 
where a,j3 = A,fi, with a cutoff at r c = \/3<Jq (<7o is 
the unit length). The choice a B = 1.2ox hampers crys- 
tallization. We impose (2a A ) 3 +2(a A + a B ) 3 + (2a B ) 3 = 
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FIGURE 1. Correlators for two translationally invariant 
quantities: the potential energy density, e, and the square of 
the density fluctuation, ^(g m j n ), at the minimal momentum al- 
lowed by the boundary conditions (1024 particles, below T mc ). 
The local swap algorithm decorrelates faster than standard MC. 

4(70. Constants C a p are chosen to ensure continuity at 
r c . The simulations are at constant volume, with parti- 
cle density fixed to o^ 3 and temperatures in the range 
[0.897r mc , 10.792T mc ], where 7/ mc is the Mode Coupling 
temperature! 13]. We use periodic boundary conditions 
on a box of size L (which discretizes momenta in units 
of q min = 2n/L) in systems with N = 512, 1024,2048 
and 4096 particles. For argon parameters, <7n = 3. 4 A, 
e/k B = 120K and r mc = 26 .4K. 

We modify the Grigera-Parisi swap algorithm f3 to 
make it local, in order to keep the algorithm in the dy- 
namic Universality Class|9] of standard Monte Carlo 
(MC). The elementary MC step is either (with proba- 
bility p) a single-particle displacement attempt or (with 
probability 1 — p) an attempt to swap particles. The swap 
is made by picking a particle at random and trying to in- 
terchange its position with that of a particle of opposite 
type, chosen at random among those at distance smaller 
than 0.6r c . Detailed balance requires that the Metropolis 
test include not only the energy variation but the change 
in the number of neighbors. The swap acceptance is inde- 
pendent of system size and grows from 0.74% at 0.9r mc 
up to 6% at 2T mc . In this work we use p = 0.5 (named 
local swap from here on) and p = 1.0 (named standard 
MC). The time unit to is N/p elementary steps. 

In fig. 

we show that the time-correlator of the lo- 
cal swap is in general different from the one of the stan- 
dard MC. The main difference is given by the absence of 
the cage effect when the particles are allowed to swap. 
This is seen by the absence of the plateau in the time 
correlator in the swap dynamics. This implies that Mode 
Coupling transition is rather ill-defined within the swap 
dynamics. In fact, at mean field level Mode Coupling sin- 
gularity is seen as due to the divergence of the caging 
time scale 1 13]. The Theory of Critical Phenomena sug- 
gests that the decay of time correlators of two different 



TABLE 1. Stretching exponent |3 of the time correlator of the energy 
at different temperatures. The fit to a stretched exponential form has been 
performed in the region C < 0. 1 
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2.or mc 


1.08 


1.02 r mc 


0.97 r mc 


0.92 r mc 


0.89 T mc 


j8 


0.99(3) 


0.99(8) 


0.98(12) 


0.89(13) 


0.83(4) 


0.93(12) 



dynamics sharing both the local nature and the conser- 
vation laws is given by the same fuction of the correla- 
tion length 4 (up to non universal numerical prefactors). 
They are said in fact to belong to the same dynamic Uni- 
versality Class 0| . Since at very long times that decay is 
expected to be exponential it is possible to define the ex- 
ponential autocorrelation time 1 15] which identifies the 
longest characteristic relaxation time t: 
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where 
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the dynamic exponent z depending only the on the Class 
of the dynamics. 

In a critical phenomenon one expects that eventually 
all the observables with the same symmetry decay with 
the same rate 1/t since they are all coupled to the slowest 
mode. The order parameter is supposed to be such slow- 
est mode and its decay is purely exponential. In fig. (0 
we show the large time decay of different translational 
invariant quantities obtained with the local swap. At high 
temperatures the exponential time is clearly the same, on 
the other hand at lower temperatures this is more diffi- 
cult to show, since the point where in log-linear repre- 
sentation all the curves become parallel moves at lower 
and lower values of C. Actually, in literature the time de- 
cay of glass-former is more often described by stretched 

exponentials like Coif) ~ exp— , with j8 < 1. The 
stretching would arise from the contribution of many re- 
gions relaxing with different time scales. For the quan- 
tities we study in this work the stretched exponential 
fit nicely the time correlators only at intermediate times 
while, as we show in tabled in the asymptotic region 
C < 0. 1 the decay is in fact purely exponential. It would 
be very difficult to get this result with standard algorithm 
since at large t the correlator is very noisy. 

The advantage of local swap over standard MC is not 
limited to the lack of the cage effect. This would be of no 
use if the exponential time of the two dynamics was the 
same. We see instead in Fig. Q that, although the T of the 
local swap seems to diverge, it remains always two order 
of magnitude smaller than the exponential time for the 
standard MC. Because of eq. the fact that this ratio is 
roughly constant at all temperatures suggests that z is the 
same for both dynamics and then that they belong to the 
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FIGURE 2. Time correlators for N = 1024 particles (N = 
2048 at T/T mc = 0.897) with the swap algorithm. Top panel: 
temperature variation of the correlator of e, close to T mc . Bot- 
tom panels: correlators of several translational invariant quan- 
tities. Asymptotically, three parallel straight lines should be 
found, with slope — 1/T exp . At the highest temperature, (bot- 
tom,left), this is clearly observed. At lower temperatures (bot- 
tom,right) we do not have enough independent measurements 
to see clearly this common slope. 
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FIGURE 3. Top: exponential time Tf as , of the translational 
invariant quantities with the swap dynamics for N = 1024 
particles (N = 2048 at T/T mc = 0.897). Bottom: the ratio of 
the exponential times for standard MC (t 5 / ow ) and local swap 
dynamics (T/ fls/ ). 



same Universality class. As anticipated in the introduc- 
tion, reducing T allows us to obtain results with an higher 
statistical accuracy. We show in fact in table [2] some pa- 
rameters describing a selection of recent numeric simu- 
lations using standard algorithms (molecular dynamics) 



TABLE 2. Statistical accuracy in the most recent 
Molecular Dynamics studies of simple glass-formers 



150000 



Paper Size* Lowest T' Time* 

Yu & Carruzzo [16] 512 1.0ir mc 10* 

Berthier|T7J 1372 0.97T mc 100 

Lacevic etal.[181 8000 1.03r mc 222 



* Number of particles. 

^ Lowest temperature where the authors equilibrated the 
system. 

** simulation length in units of relaxation time. 

* the main finding of this paper is that the specific heat 
yields a thermalization time about 100 times larger than 
the alpha correlation time. In our simulations, this is the 
relaxation time considered. 



TABLE 3. Statistical accuracy 
reached by means of the new local swap 



Temperature 


Size 


Time 


0.97r mc 


1024 


100000 


0.897r mc 


1024 


27000 


0.897r mc 


2048 


3000 


0.897r mc 


4096 


400 



to be compared in table [3] with the results that we were 
able to obtain with the local swap algorithm. 

FINITE SIZE EFFECTS 

Since the possible order parameter of the glass transition 
is unknown, a direct measure of the correlation length is 
difficult, but one may detect it indirectly through Finite- 
Size Scaling lUHl by performing numeric simulations of 
system with different sizes. In fact, in the critical region, 
quantities related with fluctuations (such as the specific 
heat) depend on the size of the system. Under the hypoth- 
esis that the order parameter belongs to the class of trans- 
lational invariant quantities, its exponential time T has to 
grow with system size as a power law (the exponent be- 
ing z) until the correlation length becomes smaller than 
the linear size of the system. Thus, in general, Eq. (|6) 
holds in the thermodynamic limit while it has to be mod- 
ified in systems of finite size fl5ll : 

T~max(L,£) z . (7) 

Another quantity of interest is the specific heat C v . As a 
matter of fact it must show a critical behavior if a sec- 
ond order phase transition is the responsible for the glass 
transition. In this case the approach to the thermody- 
namic value is rather expected to be exponential in the 
regime where L < ^ . 

In Fig. (0} we show then the relaxation time and the 
specific heat dependence on the size of the simulation 
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FIGURE 4. Top: size dependence of the exponential time at 
the two smallest temperatures studied, T/T mc = 0.921,0.897. 
Bottom: size dependence of the specific heat. 



box, L. The main result is that at the smallest studied 
temperature T =0.8977i nc , a noticeable growth of both 
quantities with L is found up to L ~ 12.7ab (~ 4.3nm 
for argon parameters) which is our estimation of the 
correlation length The quantity defined by eq. Q 
represents the largest relaxation time, which is the one 
affected by criticality. Of course, different time-scales 
(e.g. the plateau of intermediate scattering functions) 
might display a different dependence on L 11911 . 

The correlation length that we have found is com- 
parable with the experimental domain size for co- 
operative dynamics ia EH, and well above previous 
simulations! 17]. Aiming to estimate to dynamic expo- 
nent z one needs the value of | at larger temperatures, 
however up to T = 0.9217' mc no finite-size effects are de- 
tected for both quantities. We believe that one should 
study smaller system to find them, but unfortunately they 
crystallize quickly below 2.13r mc . Taking 8oo as upper 
bound for | at T = 0.921 r mc and utilizing the values of 
T at the two smallest temperature we derive from Eq. l|6) 
z~ 1.2 as a rough estimate of the upper bound for the 
dynamic exponent. 



CONCLUSIONS 

Our findings are relevant in a twofold way. 

Fast Monte Carlo algorithm We have shown that fast 
algorithms may open for numeric simulation a win- 
dow on the phenomenology of the glass transition 
that was previously available only to experiments. 
On the other hand, the local swap dynamics in- 
troduced works quite well with the simple model 
studied here but it becomes less effective when the 
parameters of the model change. In a similar soft 



sphere model |16] with Ob = \.Aga, the swap ac- 
ceptance reduces of two order of magnitudes. 
Finite size effects Having looked for a diverging corre- 
lation length we must at the moment content with a 
static correlation length that, close to the glass tran- 
sition, is an order of magnitude larger than the range 
of the interactions. This challenges the common 
view that in glasses only dynamic quantities be- 
come correlated at large scales. More importantly, 
the study of the dynamics of translationally invari- 
ant quantities appears as a challenge to experimen- 
talists. While measurements of the frequency de- 
pendence of the specific heat| 161 12211 are an appeal- 
ing possibility to estimate the potential energy re- 
laxation time, the correlation-length could be stud- 
ied by Finite-Size Scaling of the specific-heat and 
of relaxation times in films or in larger pores than 
previously used to confine glass-formers. In fact ex- 
periments in films|23] and nanopores|24] show that 
the glass transition changes in samples with one or 
more dimensions of nanometric scale 2 . Interest- 
ingly, the specific-heat of toluene confined on pores 
of diameter 8.7 nm l24ll . close to its glass tempera- 
ture, is significantly smaller than for bulk toluene, 
which could signal a correlation length well above 
the nanometric scale. 

At a qualitative level one may note that a second or- 
der phase transition is not the only scenario predict- 
ing Finite Size effects in the specific heat. In fact, 
in the random first order transition picture of the 
glass transition [25] recently revised by Biroli and 
Bouchaud 1 26] the contribution to the specific heat 
from the configurational entropy s c (T) (the loga- 
rithm of the number of metastable states) is size- 
able only when the system is large enough to con- 
tain many regions in different states. This happens 
when 

L d s c (t)>\. (8) 

To discriminate between the two interpretations a 
direct measurement of the correlation length of the 
energy would be needed. Nevertheless, in 12711 we 
have shown some evidence for a singular behavior 
of the specific heat, consistent with the second order 
phase transition interpretation. 

ACKNOWLEDGMENTS 

We thank G. Biroli for pointing out that the random first 
order transition picture may explain finite size effects in 

2 although it is difficult to disentangle Finite-Size Scaling from the 
effects of the interaction with the substrate. 



the specific heat. P. V. was supported by the EC (contract 
MCFI-2002-01262). We were partly supported by MEC 
(Spain), through contracts BFM2003-08532, FIS2004- 
05073 and FPA2004-02602. The total CPU time de- 
voted to the simulation (carried out at BIFI PC clusters) 
amounts to 10 years of 3 GHz Pentium IV. 



REFERENCES 

1. L. Chayes, V. J. Emery, S. A. Kivelson, Z. Nussinov, 
and G. Tarjus, Physica A Statistical Mechanics and its 
Applications 225,'l29-153 (1996). 

2. C. Donati, S. Franz, G. Parisi, and S. C. Glotzer, J. 
Non-Crys.Sol. 307, 215-224 (2002). 

3. H. E. Castillo, C. Chamon, L. F. Cugliandolo, J. L. Iguain, 
and M. P. Kennett, Phys. Rev. B 68, 134442-1-134442-41 

(2003) . 

4. G. Biroli, and J. -P. Bouchaud, Europhys. Lett. 67, 21-27 

(2004) . 

5. S. Whitelam, L. Berthier, and J. P. Garrahan, Phys. Rev. 
Lett. 92, 185705-1-185705-4 (2004). 

6. P. G. Debenedetti, Metastable Liquids, Princeton 
University Press, 1997. 

7. C. A. Angell, K. L. Ngai, G. B. McKenna, P. F McMillian, 
and S. W. Martin, J. Appl. Phys. 88, 3133-3156 (2000). 

8. P. G. Debenedetti, and F. H. Stillinger, Nature 410, 
259-267 (2001). 

9. J. Z. Justin, Quantum Field Theory and Critical 
Phenomena, Oxford University Press, 2002. 

10. J. Cardy, Scaling and Renormalization in Statistical 
Physics, Cambridge University Press, Cambridge, 1996. 

11. D. Amit, and V. Martin-Mayor, Field Theory, the 
Renormalization Group and Critical Phenomena, World 
Scientific Singapore, in press, 2005. 

12. L. Berthier, Phys. Rev Lett. 91, 055701-1-055701^1 
(2003). 

13. W. Gotze, and L. Sjogren, Rep. Prog. Phys. 55, 241-336 
(1992). 

14. T. S. Grigera, and G. Parisi, Phys. Rev. E 63, 045102-1- 
045102-4 (2001). 

15. A. D. Sokal, in Functional Integration: Basics and 
Applications, eds. C. DeWitt-Morette, P. Cartier and A. 
Folacci, Plenum, New York, 1997. 

16. C. C. Yu, and H. M. Carruzzo, Phys. Rev. E 69, 
051201-1-051201-10 (2004). 

17. L. Berthier, Phys. Rev. E 69, 020201-1-020201-4 (2004). 

18. N. Lacevic, F. W. Starr, T. B. SchrAyder, and S. C. 
Glotzer, J. Chem. Phys. 119, 7372-7387 (2003). 

19. K. Kim, and R. Yamamoto, Phys. Rev. E 61, R41-R44 
(2000). 

20. M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99-128 (2000). 

21. E. V. Russell, and N. E. Israeloff, Nature 408, 695-698 
(2000). 

22. N. O. Birge, and S. R. Nagel, Phys. Rev. Lett. 54, 
2674-2677 (1985). 

23. D. S. Fryer, et al., Macromolecules 34, 5627 (2001). 

24. D. Morineau, Y. Xia, and C. Alba-Simionesco, J. Chem. 
Phys. 117, 8966-8972 (2002). 

25. T. Kirkpatrick, D. Thirumalai, and P. Wolynes, Phys. Rev. 
A 40, 1045-1054 (1989). 



26. G. Biroli, and J. -P. Bouchaud, /. Chem. Phys. 121, 
7347-7354 (2004). 

27. L. A. Fernandez, V. Martin-Mayor, and P. Verrocchio, 
preprint condmat 0504327 (2005). 



